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Abstract 

Low-dimensional  flow  dynamical  systems  are  susceptible  to  instabilities  after  long¬ 
time  integration.  In  this  paper,  we  investigate  the  stability  of  such  two-dimensional 
models  constructed  from  Karhunen-Loeve  expansions  for  flows  past  a  circular  cylin¬ 
der.  We  first  demonstrate  that  although  the  short-term  dynamics  may  be  predicted 
accurately  with  only  a  handful  of  modes  retained,  instabilities  arise  after  a  few 
hundred  vortex  shedding  cycles.  We  then  propose  a  dissipative  model  based  on  a 
spectral  vanishing  viscosity  (SVV)  diffusion  convolution  operator  as  an  effective  way 
of  stabilizing  low-dinrensional  Galerkin  systems. 
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1  Introduction 


Proper  orthogonal  decomposition  (POD)  is  a  methodology  that  first  identifies 
the  few  most  energetic  modes  in  a  time-dependent  system,  and  second  provides 
a  means  of  obtaining  a  low-dinrensional  description  of  the  system’s  dynamics 
[1] .  A  particular  effective  approach  is  the  method  of  snapshots ,  first  proposed 
in  [2]  for  flow  systems,  that  makes  the  method  easy  to  implement  in  practice. 
POD  has  been  successfully  implemented  in  conjunction  with  experimental 
(e.g.,  [3-6,?])  as  well  as  with  numerical  studies  (e.g.,  [7,2,8-12])  in  thermal 
convection,  shear  layers,  cavity  flows  and  external  flows,  to  mention  just  a 
few. 
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In  some  of  the  aforementioned  studies  ad  hoc  viscosity  models  have  been  in¬ 
corporated  to  produce  stable  simulations  (e.g.,  [11]  for  cavity  flows)  while  in 
others  there  was  no  explicit  stabilization  technique  incorporated  (e.g.,  [7]  for 
cylinder  flows).  In  particular,  as  the  number  of  modes  increases  above  a  cer¬ 
tain  threshold  the  POD-based  model  seems  to  be  stable  at  least  for  short-time 
integration.  The  stability  of  the  model,  however,  is  strongly  dependent  on  the 
flow  geometry  and  regime,  while  most  of  the  fixes  for  one  flow  are  not  effective 
in  another  flow.  Also,  most  of  the  studies  so  far  have  addressed  short-term  dy¬ 
namics  while  instabilities  may  arise  after  thousands  of  convective  time  units. 

From  the  theoretical  standpoint,  it  is  well  known  that  the  system  of  ordinary 
differential  equations  derived  from  the  Galerkin  projection  of  a  dissipative 
PDE  may  be  unstable  for  the  long-term  dynamics,  see  [13].  A  potential  ap¬ 
proach  to  restore  dissipation  back  to  the  low-dimensional  system  is  nonlinear 
Galerkin  projection,  which  is  based  on  concepts  of  approximate  inertial  mani¬ 
folds,  see  ([14],  [15]  and  [16]).  For  fully  discrete  systems,  the  nonlinear  Galerkin 
method  has  been  shown  to  be  stable  for  the  long-term  dynamics  but  some  sen¬ 
sitivity  to  initial  data  was  also  revealed,  see  [17].  In  practice,  this  method  works 
effectively  as  we  have  recently  shown  in  [18]  using  a  low-dimensional  system 
constructed  from  experimental  (Particle  Image  Velocimetry)  data.  However, 
we  have  encountered  several  other  reduced  flow  model  systems  for  which  such 
stabilization  proved  inadequate.  We  will  demonstrate  this  behavior  in  the  fol¬ 
lowing. 


In  this  paper,  we  present  an  alternative  stabilization  strategy  based  on  the 
spectral  vanishing  viscosity  (SVV)  method.  SVV  was  first  introduced  in  [19] 
in  the  context  of  constructing  monotonicity  preserving  discretizations  to  hy¬ 
perbolic  conservation  laws.  More  recently,  it  has  been  employed  successfully 
in  formulating  alternative  large-eddy  simulation  (LES)  approaches  [20].  Also, 
in  [21],  the  Legendre  spectral  vanishing  method  was  shown  to  effectively  con¬ 
trol  the  Gibbs  phenomenon,  while  in  [22],  the  SVV  approach  was  employed  in 
two-dimensional  simulation  of  waves  in  stratified  atmosphere. 

The  spectral  vanishing  viscosity  approach  guarantees  an  essentially  non-oscillatory 
behavior  although  some  small  oscillations  of  bounded  amplitude  may  be  present 
in  the  solution.  This  theory  is  based  on  three  key  components: 

(1)  A  vanishing  viscosity  amplitude  which  decreases  with  the  mode  number; 

(2)  A  viscosity-free  spectrum  for  the  lower,  most  energetic  modes;  and 

(3)  An  appropriate  viscosity  kernel  for  the  high-wave  numbers. 

SVV  is  especially  suitable  for  hierarchical  discretizations  such  as  the  proper 
orthogonal  decomposition  where  global  energetically-ordered  modes  are  in¬ 
volved.  This  implies  that  SVV  preserves  the  inherent  energetic  scale  sepa¬ 
ration  while  it  also  maintains  monotonicity  of  the  total  variation  bounded 
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Fig.  1.  Computational  domain. 

(TVB)  kind  by  controlling  the  high-frequency  components.  This  effective  reg¬ 
ularization  is  determined  by  parameters  whose  range  is  given  directly  by  the 
theory  for  advection-dominated  systems.  More  recent  work  has  extended  the 
method  to  superviscosity  formulations,  first  by  Tadmor  [23]  and  later  by  Ma 
[24,25],  in  order  to  extend  the  range  of  the  viscosity-free  spectrum. 

In  the  following,  we  first  demonstrate  a  few  cases  where  instabilities  arise  and 
subsequently  we  introduce  the  SVV  method  and  present  asymptotically  stable 
results  by  incorporating  SVV. 


2  Mathematical  Formulation 


2.1  Direct  Numerical  Simulation 


We  consider  here  flow  past  a  circular  cylinder  for  which  both  two-  and  three- 
dimensional  POD  models  have  been  constructed  in  [7]  and  [26],  respectively. 
These  models  were  stable  for  tens  and  even  hundreds  of  shedding  cycles 
without  incorporating  any  stabilization  scheme.  We  will  examine  the  stabil¬ 
ity  of  these  flows;  in  particular  for  the  concepts  developed  here  we  consider 
two-dimensional  uniform  flow  past  a  circular  cylinder  at  Reynolds  number 
Re  =  100  and  Re  =  500. 

The  computational  domain  is  shown  in  figure  1.  Uniform  steady  or  time- 
dependent  boundary  conditions  are  imposed  at  the  inflow  boundary  Ti.  Uni¬ 
form  velocity  is  also  imposed  on  T3  and  T4  while  on  T2  the  zero  Neumann  con¬ 
dition  on  velocity  is  imposed.  On  the  cylinder  surface  T5  the  no-slip  boundary 
condition  is  prescribed.  Converged  solutions  were  obtained  using  the  spectral / hp 
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Fig.  2.  Instantaneous  vorticity  contours  at  Re=100  (upper)  and  Re=500  (lower). 

element  method  [27] .  Typical  results  that  show  the  differences  in  spatial  scales 
in  terms  of  vorticity  at  Re=100  and  Re  =  500  are  shown  in  figure  2. 


2.2  POD  Models 


We  employed  50  snapshots  of  DNS  data  in  order  to  construct  the  low-dimensional 
models  using  the  proper  orthogonal  decomposition.  We  briefly  review  this  pro¬ 
cedure  next. 

Let  us  decompose  the  total  flow  field  V  as 

V  (x,  t)  =  U0(x)  +  u(x,  t) 
where  U0  is  the  time-averaged  field. 

Then,  we  extract  the  POD  modes,  based  on  the  DNS  data,  which  are  eigen¬ 
vectors  of  a  covariance  matrix  C;  its  elements  are  computed  as  follows 

Cij  =  iu(x,fj) -u(x, fj)dx  =  J (u(x,y,ti)u(x,y,tj)+v(x,y,ti)v(x,y,tj)^dxdy  , 

where  u,  v  are  the  two  components  of  the  velocity  vector  u.  This  is  the  snap- 
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Fig.  3.  POD  eigenvalues  for  Re  =  100  (+)  and  Re  =  500  (*). 

shot  method  formulation ,  see  [28-30] .  The  matrix  C  represents  the  correlation 
between  temporal  points  since  the  spatial  variable  has  been  integrated  out. 
We  then  compute  the  eigenvectors  of  the  above  covariance  matrix,  denoted  by 
a,  and  the  POD  modes  denoted  by  <p(x,y).  Specifically,  the  vector  <f>(x,y )  is 
given  by 

N 

f>u(x,y)j  ^  ] Q'j(tj)'U'(x,y jtj) 

i=  1 

N 

f>v{x,y)j  =  J2aj(U)v(x,y,ti) 
i—  1 

where  N  is  the  total  number  of  snapshots,  (f)u  and  (j)v  are  the  components  of 
the  vector  (p(x,y),  and  j  is  the  mode  index.  The  corresponding  eigenvalues 
are  ordered  and  plotted  in  figure  3  for  Re  =  100  and  Re  =  500. 

We  employ  the  hierarchical  POD  modes  obtained  from  the  DNS  data  as  a  basis 
to  represent  the  velocity  held.  In  addition,  we  employ  a  Galerkin  projection 
of  the  Navier-Stokes  equations  onto  spatial  modes  to  obtain  the  system  of 
ordinary  differential  equation  that  governs  the  dynamics  of  the  system. 

We  express  the  two-dimensional  held  u  as  the  linear  combination  of  the  POD 
modes 


N 

u(x,y,t)  =  '52<f>u(x,y)jaj(t), 

3= 1 

N 

v{x,y,t )  =  Y!l<Pv(x,y)jaj(t), 

3= 1 
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where  cl  jit)  are  the  unknown  coefficients.  The  Galerkin  projection  of  the 
Navier-Stokes  equations  gives 


+  (V  •  V)V  +  Vp 


-^V2V  )  d-x  =  0 


(i) 


where  the  projection  vector  is  <fi  =  [cj)u,  (pv]T ,  extracted  from  DNS.  We  nse  the 
divergence-free  eigenmodes  so  the  pressure  term  inside  the  domain  is  elimi¬ 
nated  via  integration  by  parts.  Also  the  Dirichlet  and  the  outflow  conditions 
imposed  at  the  boundaries  lead  to  the  vanishing  of  contributions  from  the 
pressure  on  those  boundaries  in  the  integration  by  parts  procedure. 


The  Galerkin  projection  leads  to  the  dynamical  system: 


8 CLj  (t) 

dt 


/(a) 


(2) 


with  a  =  [ai,  a-2,  ■  ■  ■]■  The  term  /(a)  includes  the  convective  and  viscous  terms 
and  has  the  form: 


/(a)  =  ~(J  ■  (0i0fc)dx)  aittk 

~  (Jfe  )  V(pJV(f)ldx  +  j  <f>.V  ■  (0iUo)dx  +  j  ■  (Uo^Jdx^  a* 
-  ('/'  0,-V  •  (UoUo)dx  +  ^  j  V^VUqcZx)  . 

In  the  following  we  investigate  the  time  evolution  of  the  modal  coefficients 
j  =  1,2,...;  we  also  refer  to  aj  as  the  “mode”  j. 


3  Instability  of  POD  Flow  Models 


First  we  present  results  from  the  long-time  integration  of  the  Re  =  100  case 
with  steady  uniform  inflow.  It  was  found  in  [7]  that  a  6-mode  POD  system 
gives  accurate  results  in  comparisons  with  the  original  DNS  data,  at  least 
for  the  short-time  dynamics.  Indeed,  various  reduced  models  we  tested  again 
with  N  >  6  are  stable  after  short-time  integration,  and  in  fact  for  times  up  to 
several  hundreds  of  shedding  cycles.  Our  experiments,  however,  show  that  all 
models  are  asymptotically  unstable.  For  example,  a  6-mode  model  is  steady 
for  to  to  40  shedding  cycles  (about  200  convective  time  units)  as  shown  in 
figure  4  but  it  diverges  for  time  t  >  200.  As  the  number  of  modes  increases 
the  stability  of  the  model  is  enhanced.  So,  a  10-mode  model  is  stable  for  up  to 
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6  mode 


6  mode 


Fig.  4.  Re=100:  Time  history  showing  the  onset  of  instability  (bifurcation)  in  the 
6-mode  POD  model. 


500  shedding  cycles,  i.e.  more  than  3,000  convective  time  units,  but  eventually 
all  modes  diverge  as  shown  in  figure  5.  In  particular,  a  bifurcation  to  a  new 
asymptotic  state  is  shown.  However,  for  the  first  50  shedding  cycles  the  POD 
predictions  are  in  very  good  agreement  with  the  DNS  data  as  shown  in  figure 
6. 


The  exact  onset  of  instability  depends  on  the  number  N  of  modes  retained  in 
the  reduced  model  as  well  as  on  the  Reynolds  number.  At  the  higher  Reynolds 
number  (Re  =  500)  the  onset  of  instability  arises  earlier  even  for  a  higher- 
order  model.  For  example,  in  figure  7  we  show  the  time  history  of  the  modes 
for  a  20-mode  POD  model.  The  instability  here  sets  in  at  about  100  shedding 
cycles  into  the  time  integration.  This  result  is  typical  of  several  other  models 
we  constructed  for  the  uniform  steady  inflow. 

However,  not  all  low-dimensional  systems  are  asymptotically  unstable.  Our 
experiments  show  that  forced  systems ,  i.e.  systems  with  an  imposed  time  scale 
through  external  forcing,  may  be  stable  at  all  times.  To  this  end,  we  consider 
the  flow  past  a  cylinder  again  at  Re  =  500  but  with  a  small  sinusoidal  velocity 
component  added  at  the  inflow,  with  10%  amplitude  forced  at  the  Strouhal 
frequency.  The  resulting  POD  system  predicts  the  expected  lockin  state,  in 
agreement  with  DNS,  and  it  is  stable  [31].  In  figure  8  we  show  the  time  history 
for  the  same  set  up  and  parameters  as  in  figure  7.  It  shows  stability  at  the 
time  of  the  instability  onset  of  the  flow  described  in  7  but  also  at  much  longer 
times  (not  shown  here).  The  Galerkin  model  in  the  unsteady  inflow  case  is 
based  on  a  modification  of  the  system  of  equations  (2)  to  include  a  penalty 
term  that  facilitates  the  time-dependent  boundary  conditions  in  the  reduced 
POD  system. 
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Fig.  5.  Re=100:  Time  history  showing  the  onset  of  instability  (bifurcation)  in  the 
10-mode  POD  model. 

4  The  Spectral  Vanishing  Viscosity  Model 


Tadmor  [19]  first  introduced  the  concept  of  spectral  vanishing  viscosity  (SVV) 
using  the  inviscid  Burgers’  equation 


d_ 

dt 


u(x,t )  + 


d 

dx 


u2(x,  t) 
2 


0, 


(3) 


subject  to  given  initial  and  boundary  conditions.  The  distinct  feature  of  solu¬ 
tions  to  this  problem  is  that  spontaneous  jump  discontinuities  (shock  waves) 
may  be  developed,  and  hence  a  class  of  weak  solutions  can  be  admitted.  Within 
this  class,  there  are  many  possible  solutions,  and  in  order  to  single  out  the 
physically  relevant  one  an  additional  entropy  condition  is  applied,  of  the  form 

d_  f u2(x,t)\  d( u3(x,t)\ 
dt\  2  j  3^  3  )  ~ 


In  low-dimensional  systems  it  has  been  found  that  unstable  behavior  is  asoci- 
ated  with  multiple  spurious  steady  states  [13],  and  this  is  consistent  with  the 


10  mode 


10  mode 


Fig.  6.  Re=100:  Phase  portrait  for  the  first  50  shedding  cycles  for  the  10-mode 
model.  The  lines  denote  POD  predictions  and  the  triangles  DNS  data. 


above  observation.  Tadmor  (1989)  introduced  the  spectral  vanishing  viscos¬ 
ity  method,  which  adds  a  small  amount  of  mode- dependent  dissipation  that 
satisfies  the  entropy  condition,  yet  retains  spectral  accuracy.  It  is  based  on 
viscosity  solutions  of  nonlinear  Hamilton- Jacobi  equations,  which  have  been 
studied  systematically  in  [32],  Specifically,  the  viscosity  solution  for  the  Burg¬ 
ers’  equation  has  the  form 


d 

dt 


u(x,t )  + 


d_ 

dx 


u2(x,  t) 
2 


8  <9w 

6  dx  6  dx 


(5) 


where  e(—>  0)  is  a  viscosity  amplitude  and  Qe  is  a  viscosity  kernel.  Convergence 
may  then  be  established  by  compactness  estimates  combined  with  entropy 
dissipation  arguments  [19].  To  respect  spectral  accuracy,  the  SVV  method 
makes  use  of  viscous  regularization  and  equation  (5)  may  be  rewritten  in 
discrete  form  (retaining  N  modes)  as  in  our  POD  model 


d  d 

mUNM  +  ai 


>N  (VQM)y 


Qn  * 


duN 

dx 


(6) 
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Fig.  7.  Re=500;  steady  inflow :  Time  history  showing  the  onset  of  instability  (bifur¬ 
cation)  in  the  20-mode  POD  model. 

where  the  star  (*)  denotes  convolution  and  Vn  is  a  projection  operator.  QN  is 
a  viscosity  kernel,  which  is  only  activated  for  high  wave  numbers.  In  Fourier 
space,  this  kind  of  spectral  viscosity  can  be  efficiently  implemented  as  multi¬ 
plication  of  the  Fourier  coefficients  of  un  with  the  Fourier  coefficients  of  the 
kernel  QN,  i.e., 


d 


e 


dx 


Qn  * 


Dun 

dx 


-e  Y,  k2Qk(t)uk(t)elkx, 

M<\k\<N 


where  k  is  the  wave  number,  N  the  number  of  Fourier  modes,  and  M  the 
wavenumber  above  which  the  spectral  vanishing  viscosity  is  activated.  In  the 
POD  context,  we  also  assume  that  this  implementation  of  convolution  is  valid 
in  the  modal  space. 


Originally,  Tadmor  (1989)  used 


j  0,  |  k  \<  M 
[l,\k\>M, 


(7) 
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Fig.  8.  Re=500;  oscillatory  inflow :  Time  history  showing  stability  for  the  20-mode 
POD  model. 


with  eM  ~  0.25  based  on  the  consideration  of  minimizing  the  total- variation 
of  the  numerical  solution.  In  subsequent  work,  however,  a  smooth  kernel  was 
used,  since  it  was  found  that  the  C°°  smoothness  of  Qk  improves  the  resolution 
of  the  SVV  method.  For  Legendre  pseudo-spectral  methods,  Maday  et  al.  [33] 
used  e  ~  iV-1,  activated  for  modes  k  >  M  rs  5 y/N,  with 

„  (fc-JV)2 

Qk  =  e  C*-M)a  t  k  >  M.  (8) 


In  order  to  see  the  difference  between  the  convolution  operator  on  the  right- 
hand-side  in  equation  (6)  and  the  usual  viscosity  regularization,  following  Tad- 
mor  [34],  we  expand  as 


d_ 

dx 


Qn  * 


8un 

dx 


82un 
e  dx2 


[Rn(x,  t) 


* 


duN 

dx 


(9) 
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Fig.  9.  Normalized  viscosity  kernels  for  the  spectral  vanishing  viscosity  (dash  line 
(7  =  0  and  solid  line  (7  =  5)  and  the  Kraichnan/Chollet-Lesieur  viscosity  (dash-dot 
line) . 

where 


N 

RN{x,t)=  £  Rk(t)eikx;  Rk(t ) 

k=—N 


1  -Qk{t)  \k\>M 
1  \k\  <  M 


(10) 


The  extra  term  appearing  in  addition  to  the  first  standard  viscosity  term 
makes  this  method  different.  It  measures  the  distance  between  the  spectral 
(vanishing)  viscosity  and  the  standard  viscosity.  This  term  is  bounded  in  the 
L2  norm  similarly  to  the  spectral  projection  error.  In  this  paper  we  refer  to 
the  viscosity  as  vanishing  as  the  theory  requires  that 


1 

NelogN ’ 


0  <  1 


and  thus  e  — »  0  for  the  high- resolution  limit. 


At  this  point  it  is  also  instructive  to  compare  the  spectral  vanishing  viscosity 
to  the  aforementioned  spectral  eddy-viscosity  introduced  by  Kraichnan  [35] 
and  modified  by  Chollct-Lesieur  [36,37].  The  latter  has  the  non-dimensional 
form  [37] 

v{k/N)  =  /1q"3/2  [0.441  +  15.2exp(—3.03N/  k)\,  K0  =  2.1  (11) 


Comparing  the  Fourier  analog  of  this  eddy-viscosity  employed  in  LES  [36] 
to  the  viscosity  kernel  Qk(k,  M,  N)  introduced  in  the  SVV  method,  figure  9 
shows  both  viscosity  kernels  normalized  by  their  maximum  value  at  k  =  N . 
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For  SVV  two  different  values  of  the  cut-off  wavenumber  are  considered,  i.e., 

M  =  CVN  for  C  =  0  and  C  =  5,  (12) 

and  are  shown  in  the  plot  of  figure  9.  In  particular,  the  solid  line  can  be 
thought  of  as  a  stability  barrier  above  which  monotonicity  and  thus  stability 
is  not  guaranteed.  On  the  other  hand,  the  dash  line  can  be  thought  of  as 
an  accuracy  barrier  below  which  the  convergence  of  the  method  is  affected. 
This  range  has  been  used  in  most  of  the  numerical  experiments  so  far  (see 
for  example  [33], [22], [20])  and  is  consistent  with  the  theoretical  results  [19].  In 
the  plot  it  is  shown  that,  in  general,  the  two  forms  of  viscosity  have  similar 
distributions  but  the  SVV  form  does  not  affect  the  first  one-third  or  one- 
half  of  the  spectrum  (viscosity-free  portion)  and  it  increases  faster  than  the 
Kraichnan/Chollet-Lesieur  eddy- viscosity  in  the  higher  wave  numbers  range, 
e.g.  in  the  second-half  of  the  spectrum. 


The  implementation  of  the  SVV  in  the  POD  models  (equation  (2))  is  similar 
to  the  implementation  of  Fourier  methods  presented  above  or  the  spectral/hp 
element  discretization  in  [20].  In  particular,  the  system  of  ordinary  differential 
equations  is  enhanced  as  follows 


^  _  Ma),  (13) 

where  /(a)  has  the  form  presented  in  equation  (3),  and  h( a)  contains  the 
viscosity  convolution  kernel,  i.e. 


h(  a)  =  eQj 


<9Uq  d(P:) 
dx  dx 


N 

dx  +  y~)  a,i(t) 

i= 1 


dx  dx 


(14) 


In  this  derivation,  integration  by  parts  is  used  and  the  fact  that  boundary 
contributions  vanish  because  of  the  specific  boundary  conditions  employed. 
In  view  of  equation  (14),  we  can  see  that  only  the  higher  modes,  i.e.  mode 
numbers  greater  than  M,  will  be  affected  by  the  viscosity  kernel. 


In  order  to  contrast  the  effect  of  SVV  to  simply  adding  artificial  viscosity  in 
the  POD  model  we  have  also  performed  simulations  with  model;  the  total 
(non-dimensional)  viscosity  in  this  case  is  mode-dependent  and  given  by 


Vk  = 


Re 


k  <  M 


1  (  k  —  N  \2 

+  Oexp  O-aP  k  >  M 


(15) 


where,  M  is  the  cut-off  wavenumber  as  before  and  C  is  a  constant. 
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Re(modes  used) 

Artificial  Viscosity 

Spectral  Vanishing  Viscosity 

Cut-off  (M) 

Constant  (C) 

Cut-off  (M) 

Constant  (e) 

100  (6  modes) 

4 

0.008333 

4 

0.0137 

100  (10  modes) 

6 

0.04 

6 

0.00822 

500  (6  modes) 

2 

0.083333 

2 

0.24167 

500  (12  modes) 

8 

0.166666 

10 

0.12954961443 

500  (20  modes) 

n/a 

n/a 

16 

0.0320520703 

Table  1 

Parameters  used  in  the  stabilization  schemes. 

This  artificial  viscosity  model  is  chosen  to  have  the  exponential  effect  that 
SVV  has,  so  to  that  the  higher  modes  affected  the  most. 


5  Stability  of  SVV-POD  Flow  Models 


Here  we  demonstrate  the  stabilization  effect  of  SVV  by  revisiting  the  flow 
examples  already  presented  in  section  2  for  cylinder  flow.  The  parameters  of 
the  SVV  model  were  chosen  guided  by  the  theoretical  estimates  and  also  by 
obtaining  the  best  agreement  with  the  original  data  for  the  first  50  shedding 
cycles.  A  summary  of  the  parameters  used  is  presented  in  table  1. 

We  first  consider  the  case  Re  =  100  and  examine  the  phase  portraits  ob¬ 
tained  for  all  modes  for  the  first  1000  shedding  cycles.  The  model  with  the 
smallest  number  of  modes  that  predicts  accurately  the  short-term  dynamics 
corresponds  to  N  =  6.  However,  over  the  span  of  1000  shedding  cycles  there 
is  loss  of  stability  of  this  specific  model  as  mentioned  earlier  (see  figure  4).  In 
figure  10  we  show  the  phase  portrait  corresponding  to  the  first  1000  shedding 
cycles;  clearly  the  initial  good  agreement  with  the  DNS  data  is  eventually 
lost.  Stability  is  recovered,  however,  if  the  SVV  correction  is  incorporated.  In 
figure  11  we  plot  again  the  phase  portrait  (over  1000  shedding  cycles)  of  the 
6-mode  system  with  the  SVV  model  included.  We  see  that  a  stable  limit  cycle 
is  established  in  agreement  with  the  DNS  data. 

Similar  results  are  valid  for  the  models  with  higher  modes.  Here,  we  compare 
the  two  viscosity  formulations  to  illustrate  what  is  different  with  SVV.  In  figure 
12  we  first  show  the  results  for  the  10-mode  POD  system  using  the  artificial 
viscosity  for  stabilization.  In  figure  13,  we  present  the  same  results  but  with 
the  SVV  model  employed.  We  see  that  the  agreement  of  the  POD  predictions 
with  the  DNS  original  data  is  uniformly  good  for  all  modes  in  contrast  with 
the  artificial  viscosity  model.  The  latter  shows  a  small  divergence  which  is 
even  greater  at  longer  times  (not  shown  here).  Also,  the  higher  modes  are 
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6  mode 


6  mode 


Fig.  10.  Re=100:  Phase  portrait  for  the  6-mode  model  without  explicit  dissipation. 
The  lines  denote  POD  predictions  and  the  triangles  DNS  data. 


6  mode 


6  mode 


Fig.  11.  Re=100:  Phase  portrait  for  the  6-mode  model  with  SVV.  The  lines  denote 
POD  predictions  and  the  triangles  DNS  data. 
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Fig.  12.  Re=100:  Phase  portrait  for  the  10-mode  model  using  the  artificial  viscosity 
model.  The  lines  denote  POD  predictions  and  the  triangles  DNS  data. 

excessively  damped.  In  order  to  compare  the  stabilization  of  this  technique 
with  the  stabilization  effect  of  the  nonlinear  Galerkin  method,  we  plot  in  figure 
14  corresponding  results  from  a  nonlinear  Galerkin  model,  see  [38]  for  details. 
The  nonlinear  Galerkin  model  is  constructed  based  on  six  dominus  (also  known 
as  masters)  and  four  servus  (also  known  as  slaves)  POD  modes.  We  see  from 
the  phase  portrait  that  this  approach  is  inadequate  in  obtaining  asymptotically 
stable  results  although  it  has  improved  the  results  compared  to  Galerkin-only 
projection.  The  same  conclusions  have  been  obtained  from  results  with  several 
other  dominus-servus  combinations  not  shown  here. 


For  the  higher  Reynold  number,  Re  =  500,  for  which  the  POD  Galerkin 
system  bifurcates  at  earlier  times  SVV  can  effectively  stabilize  the  simulation. 
In  figure  15  we  plot  the  phase  portrait  of  the  first  nine  modes  for  the  first 
1000  shedding  cycles  in  the  simulation.  We  see  that  a  limit  cycle  is  predicted 
in  excellent  agreement  with  the  DNS  data.  To  appreciate  the  effect  of  SVV  in 
the  current  simulation  we  also  present  in  figure  16  the  corresponding  results 
without  SVV  from  the  pure  POD  Galerkin  system  for  the  same  time  period, 
which  clearly  diverges.  In  this  case  the  cut-off  mode  was  set  to  M  =  16 
and  the  viscosity  kernel  e  =  C/N  was  set  as  shown  in  table  5.  Using  the 
artificial  viscosity  approach  to  stabilize  the  POD  Galerkin  model  does  not 
lead  to  accurate  results. 
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Fig.  13.  Re=100:  Phase  portrait  for  the  10-mode  model  using  the  SVV  model.  The 
lines  denote  POD  predictions  and  the  triangles  DNS  data. 

In  order  to  also  investigate  the  effect  of  SVV  on  the  higher  modes  we  plot 
separately  in  figures  17  and  18  the  phase  portraits  up  to  the  17-th  and  20-th 
mode,  respectively.  We  see  that  after  the  cut-off  mode  M  —  16  some  inaccu¬ 
racies  are  introduced,  which  are  more  pronounced  in  the  modes  18,  19  and  20. 
However,  the  amplitude  of  those  modes  is  bounded  at  all  times  in  contrast  to 
the  high  POD  modes  of  the  reduced  system  without  stabilization. 

Finally,  we  note  that  even  lower  dimensional  systems  at  Re  =  500  with  trun¬ 
cations  corresponding  to  N  =  6  and  N  —  12  are  stable  and  give  accurate 
results  for  the  SVV  parameters  shown  in  table  5. 


6  Summary  and  Discussion 


We  have  developed  a  new  approach  to  stabilizing  reduced  order  models  de¬ 
rived  from  Galerkin  projections  of  evolution  equations.  Specifically,  here  we 
have  considered  the  external  flow  past  a  cylinder  and  investigated  the  stability 
of  the  limit  (shedding)  cycle  obtained  from  a  POD-based  Galerkin  system  at 
two  values  of  Reynolds  number.  We  have  found  that  in  the  current  long-term 
time  integration  employed  all  Galerkin  models  are  asymptotically  unstable. 
However,  in  previous  similar  studies,  where  short  to  modest  length  time  in- 
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Fig.  14.  Re=100:  Phase  portrait  for  the  10-mode  model  using  the  nonlinear  Galerkin 
method.  The  lines  denote  POD  predictions  and  the  triangles  DNS  data. 

tegration  was  involved,  it  was  assumed  that  such  models  are  asymptotically 
stable.  This  instability  does  not  always  manifest  itself  as  an  explosive  growth 
that  leads  to  blow-up  but  as  a  bifurcation  that  leads  to  another  limit  cycle. 
The  precise  onset  of  this  bifurcation  depends  on  the  number  of  modes  retained 
in  the  model  and  the  Reynolds  number  as  well  as  the  flow  geometry. 

The  classical  way  of  incorporating  artificial  viscosity  to  stabilize  low-order 
models  does  not  guarantee  stability  and  affects  greatly  the  accuracy  of  the  so¬ 
lution.  On  the  other  hand,  nonlinear  Galerkin  projection,  although  potentially 
effective  for  the  right  combination  of  dominus-servus  modes,  is  not  robust;  in 
the  cases  we  considered  here  it  simply  prolonged  the  onset  of  the  instability. 
The  spectral  vanishing  viscosity  (SVV)  method  we  implemented  is  an  effec¬ 
tive  and  robust  stabilization  scheme,  ft  employs  a  convolution  viscosity  kernel, 
which  is  parametrized  by  a  viscosity  amplitude  and  a  cut-off  wavenumber. 
The  theory  does  not  give  the  precise  values  of  these  parameters  but  provides 
a  stability  range  in  terms  of  barriers.  Above  a  certain  value  of  the  cut-off  no 
stability  is  guaranteed  whereas  below  a  low  threshold  the  accuracy  of  the  low 
most  energetic  modes  is  affected.  The  viscosity  amplitude  scales  inversely  pro¬ 
portional  to  the  mode  number,  i.e.  oc  C/N ,  where  the  constant  C  is  problem 
dependent.  Here  it  was  chosen  on  the  basis  of  matching  the  model’s  short-term 
dynamics  with  the  original  data. 
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Fig.  15.  Re  =  500:  Phase  portrait  of  the  first  9  modes  in  the  20-mode  SVV-POD 
system  (lines)  compared  to  the  DNS  data  (symbols). 


We  have  defined  here  asymptotic  stability  as  stability  for  the  first  1000  shed¬ 
ding  cycles  somewhat  arbitrarily.  It  is  possible  that  after  a  much  longer  time 
interval  an  instability  is  developed.  To  this  end,  we  investigated  this  question 
and  found  that  for  the  6-mode  system,  which  can  be  integrated  easily  for  very 
long  time,  stability  was  achieved  for  more  than  1.25  million  convective  time 
units.  For  the  other  models  some  small  divergence  was  detected  after  sev¬ 
eral  thousands  of  shedding  cycles  but  a  small  change  in  the  SVV  parameters 
could  improve  the  results.  Clearly,  more  theoretical  work  is  needed  towards 
this  direction  to  provide  simple  guidelines  to  the  practitioners  of  this  method. 

The  SVV  nonlinear  stability  theory  is  based  on  the  treatment  of  the  inviscid 
Burgers  equation  originally  proposed  by  Tadmor  [19]  for  Fourier  discretiza¬ 
tion.  We  can  justify  its  use  in  the  current  context  only  heuristically  and  have 
been  motivated  by  success  in  other  applications  [22,20].  However,  a  rigorous 
justification  for  low- dimensional  models  derived  from  Galerkin  projections  is 
currently  missing  and  thus  we  do  not  have  much  insight  into  the  effectiveness 
of  SVV.  The  extra  term  appearing  in  equation  (9),  in  addition  to  the  standard 
viscosity  term,  is  perhaps  the  key  but  its  optimum  form  may  depend  on  the 
specific  dissipative  PDE  considered.  Future  work  should  address  these  issues, 
and  also  investigate  the  dissipation  spectrum  in  detail  for  larger  systems  with 
higher  number  of  modes  so  that  a  sufficient  dissipation  range  exists. 
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Fig.  16.  Re  =  500:  Phase  portrait  of  the  first  9  modes  in  the  20-mode  POD  only 
system  (lines)  compared  to  the  DNS  data  (symbols). 
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